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Abstract 

We consider a two-layer Heisenberg antiferromagnet which can be either in 
the Neel-ordered or in the disordered phase at T = 0, depending on the ratio 
of the intra- and interlayer exchange constants. We reduce the problem to an 
interacting Bose-gas and study the sublattice magnetization and the trans- 
verse susceptibility in the ordered phase, and the spectrum of quasiparticle 
excitations in both phases. We compare the results with the spin-wave theory 
and argue that the longitudinal spin fluctuations, which are not included in 
the spin-wave description, are small at vanishing coupling between the layers, 
but increase as the system approaches the transition point. We also compute 
the uniform susceptibility at the critical point to order 0(T 2 ), and show that 
the corrections to scaling are numerically small, and the linear behavior of Xu 
extends to high temperatures. This is consistent with the results of the recent 
Monte-Carlo simulations by Sandvik and Scalapino. 
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I. INTRODUCTION 



In the past few years, there has been a significant interest in the physics of quantum 
phase transitions in 2D spin systems JT]-0]. The purpose of the present communication is to 
study in detail the disordering transition in a two-layer S = 1/2 Heisenberg antiferromagnet 
described by 

H — J\ ^ S at iS a j + J2 ^ S\ y iS2,i ■ (1) 
<i,j>,a i 

Here a = 1,2; the first sum runs over nearest-neighbors, and the exchange couplings are 
assumed to be positive (see Fig. [[]). 

For small J2/J1, the model describes two weakly interacting 2D Heisenberg antiferro- 
magnets. Each of them is ordered at T = and possesses Goldstone excitations related to 
a spontaneous breakdown of a rotational symmetry. In the opposite limit, J2/J1 ^> 1, pairs 
of adjacent spins from different layers form spin singlets separated from triplet states by a 
gap, ~ J2. The presence of a gap implies that the rotational symmetry is not broken. Thus 
one should expect a disordering phase transition at some critical ratio of J2I J\- 

The two-layer Heisenberg model has attracted a lot of interest in the last few years [^-[13 



This interest was stimulated in part by the experimental observation that some of the high- 
T c superconductors contain pairs of C11O2 layers which are separated from other layers by 
a charge reservoir In addition, a two-layer antiferromagnet is probably the simplest 

non frustrated spin system which displays a quantum disordering transition of the 0(3) 
universality class. Several quantitative predictions about the behavior of observables near 
such a transition have been made recently 0, and a two-layer antiferromagnet is an ideal 
candidate to test these predictions. 

The phase diagram of eq. ([!]) has been studied numerically, by Quantum Monte-Carlo (12 



and series expansion [|T(J] techniques, and analytically, using spin-wave and mean-field 
Schwinger-boson theory There are several issues which emerged from these studies. 

Some of them are related to the universal ratios of various observables and are discussed 
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elsewhere ||14j| . Here we will focus on the properties of the system at T = 0, and on 
the corrections to scaling at finite T. The key issue we want to address at T = is 
the applicability of perturbative and self-consistent spin-wave approaches (the latter is very 
similar to the Schwinger-boson mean-field theory). It is well known that spin- wave expansion 
works extremely well for a single-layer S = 1/2 antiferromagnet. At the same time, for a two- 
layer system, spin-wave and Schwinger-boson theories yield results which are inconsistent 
with numerical simulations. In particular, the Schwinger-boson mean-field theory yields a 
critical value of the interlayer coupling ( J2/ Ji) C r ~ 4.5 |Tl| , which is nearly two times larger 
than (J2I 'Ji)cr ~ 2.5 obtained in series expansion Jl0| (J| r = 2.56) and Quantum Monte- 
Carlo [0 ( = 2.51) studies. A self-consistent spin-wave theory (see Sec [I A|) also predicts 
a very large value of (J 2 / ' J\) C r ~ 4.3. The spin-wave velocity at the critical point, c sw ~ 2J X 
is also inconsistent with the Monte-Carlo data which yield c sw ~ 1.7 — I.8J1 (where we set 
the lattice constant ao equal to unity). We will argue that the discrepancies between the 
spin-wave results and the numerical simulations have a physical origin and are related to 
the fact that in the spin-wave approach one neglects longitudinal spin fluctuations. Our 
analytical approach to the problem is based on the introduction of a triplet of S = 1 bosons 
for a pair of S = 1/2 spins (see eq. (|7|) below). In the disordered phase, this triplet of bosons 
describes the excitations above the singlet ground state of a pair, while in the ordered phase, 
where we introduce a condensate for one type of boson, the excitations are split into two 
transverse and one longitudinal magnon modes. We will show that the contributions from 
longitudinal fluctuations to J| r and c sw are substantial, which makes the 1/S expansion 
inapplicable. However, we will also show that as J 2 decreases, the spin-wave approximation 
becomes more and more reliable, and at vanishing J 2 , longitudinal spin fluctuations do not 
contribute to the sublattice magnetization and susceptibility. 

Furthermore, we will discuss the temperature dependence of the uniform susceptibility, 
Xui at the transition point. Monte-Carlo simulations have shown that the universal, linear 
behavior of Xu at J 2 = J 2 extends up to very high T ~ J\. For comparison, in a single-layer 



antiferromagnet, the deviations from linearity become substantial already at T ~ 0.6 J ||15|| . 



To understand this result, we will compute the leading nonuniversal, 0(T 2 ) correction to 
the susceptibility and show that it is numerically quite small for all T < J\. 

We start in the next subsection with the spin- wave calculations for eq. ([]]). In Sec [IT], we 
will introduce the transformation to bosons and consider in a systematic way the excitations 
in the disordered phase, the critical value of J2, and the spin- wave velocity at the critical 
point. In Sec we extend the approach to the ordered state by introducing a condensate 
for one of the bosonic fields. We will show how the triplet of excitations splits into two 
gapless transverse modes and a longitudinal mode with a finite gap. We will obtain the 
T = sublattice magnetization and the uniform susceptibility at arbitrary J2 and show how 
they deviate from spin- wave results for increasing J 2 . Finally, in Sec |IV|, we will compute 
the uniform susceptibility, Xu(T) at the critical point and show that the lattice-dependent, 
0(T 2 ) corrections to the scaling form of Xu remain small up to T = J\. Our conclusions are 
stated in Sec |V|. 

A. Spin-wave calculations 

We start our considerations with a brief discussion of the spin-wave calculations. At small 
J2, the spins are ordered antiferromagnetically in the layers and also between the layers. 
Introducing bosons via the Holstein-Primakoff transformation, and performing standard 
manipulations, we obtain two branches of spin- wave excitations with the dispersion ei(k) = 
€2(k + 7r) = efc, where to order 1/S, 

e k = U X S[{1 - vl) + (J 2 /2Ji)(l - ^)] 1/2 , (2) 

and Uf. = (cosk x + cosk y )/2. It is not difficult to show that the fluctuations near k = n are 
in-phase fluctuations of the spins in the two layers, while those near k = correspond to 
out-of-phase fluctuations. There is indeed a Goldstone mode in €fc at k — (-7T, 7i") because of 
a spontaneous symmetry breaking. The renormalized J\ and J2 differ from the couplings in 
([!]) due to the 1/S* corrections: 
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where 



1 4^ + J 2 S 1 1 (4Jxg^ + J 2 SK 1 v A^Suk + J 2 g 

(4) 

The summation in @ is over the whole Brillouin zone. The sublattice magnetization to order 
1/Sis given by N = S—5i . Evaluating 8\ with bare couplings Ji j2 , as it is required in the 1/S 
expansion, we obtain that 5± reaches a value of S = 1/2 only at a very large J 2 /Ji ~ 13.6. 
A somewhat better, though less justified estimate of J 2 r can be obtained if one formally 
considers the expressions for the renormalized couplings as self-consistent equations, and 
solve them for S = 1/2. This procedure is similar to the mean-field Schwinger-boson theory, 
and the results we obtained are similar to those of Millis and Monien 111 111: the sublattice 



magnetization first increases with J 2 , passes through a maximum, and then decreases (see 
Fig. D. There is a weak first-order disordering transition at J% m 4.36Ji. Still, the critical 
J 2 is much larger than J| r w 2.5 Ji obtained in numerical simulations. We also computed the 
spin- wave velocity to second order in 1/S, and obtained after straightforward but somewhat 



tedious calculations [18| 



^=2^2^1 + ^(1 + ^), (5) 

where now Ji j2 are the solutions of (§) and ([|) to order 1/S* 2 , and Q is a cumbersome function 
of J 2 /Ji whose explicit form we do not present. At J 2 = 0, we obtained Q ~ 0.022 which 
completely agrees with the results of other studies [ 2(J . The spin- wave velocity remains finite 



at the critical point, and it is therefore reasonable to compute it at (J 2 /Ji) C r suggested by 
numerical simulations. Assuming for definiteness that J| r = 2.55Ji, we obtained Q ~ 0.044. 
Evaluating then Ji j2 and substituting them into (§), we find c sw ~ 3.62 JiS(l + 0.094/2S 1 + 
0.026/(2S) 2 ). Observe that the 1/S 2 correction is very small. For S = 1/2, we obtain 
c sw ~ 2.03 Ji. As we mentioned earlier, this value is somewhat larger than c sw ~ 1.7— I.8J1 



extracted from the fit of the Monte-Carlo data for the uniform susceptibility |12| to the 
scaling formula. 

The main weakness of the spin-wave theory is that it assumes that the long-range or- 
der is well established, and only includes transverse spin fluctuations. However, at the 
critical point, transverse and longitudinal fluctuations become indistinguishable and should 
be treated on equal ground. We therefore proceed now to perturbative calculations which 
explicitly take the longitudinal spin fluctuations into account. 

II. DISORDERED PHASE 

The key starting point of our consideration is an observation that for sufficiently large 
J2, pairs of adjacent spins from the two planes form spin singlets. The excited state of a 
given pair is a three-fold degenerate triplet state. It is then natural to introduce a triplet of 
bosons for any given pair. Each boson describes the transformation from a singlet state to 
one of the states with S — 1. Specifically, we introduce 

Mi=Si : i + Sz,i, Li = Si t i — 5 , 2,i, (6) 

and three bosonic fields as 

Ml = af a, - b+k, LI = - (c+Ui + U iCi ), 

M+ = V2(aia - cfh), Lf = V2(atUi + UA), 

M~ = V2(cfa,i — bfci), L~ = V2(btUi + U iai ). (7) 



where C/j = y 1 — aJa, L — b^bi — cJCi. It is easy to check that the commutation relations 
for M and L are the same as for a vector and a generator of rotations: [M a , M^\ = 
ieapyM 1 ; [L a , L®\ = ie a p 1 M' 1 ; [M a ,L^] = ie^L 1 '. This in turn implies that the spin 
commutation relations for Si and S2 are satisfied. The U term, however, imposes the con- 
straint that only one boson can be excited at each lattice site. This indeed follows from the 
fact that there are only four physical states for a given pair of spins. For the physical states, 
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we have Sf = 3/4 as it should be. Notice that a similar restriction on the number of bosons 
holds also for the conventional Holstein-Primakoff transformation for S = 1/2. In this sense, 
the transformation above can be viewed as an extension of the Holstein-Primakoff transfor- 
mation to nonmagnetic states. One can also introduce an analog to the Dyson- Maleev 
transformation, but we found that the latter is less convenient for practical purposes. 

Furthermore, a conventional way to perform spin-wave calculations for a Neel-ordered 
state of a S = 1/2 system is to extend a model to large S, perform 1/S expansion, and 
set S = 1/2 at the very end of the calculations. We will now do the same for a disordered 
state. To this end, we modify the transformation to bosons by introducing a factor A <C 1 



into the square root as = y 1 — X(afa,i + bfbi + cfc/), and simultaneously introducing 
an overall factor 1\/X into all three components of Lj. It is not difficult to check that the 
commutation relations between L and M (and, hence, the spin algebra) do not change under 
this transformation; however, the value of the spin on each site in the ground state is now 
0(1/ X) ^> 1. Below, we perform a systematic perturbative expansion in A. The physical 
results, however, correspond only to A = 1. 

Eq. (|7|) has been applied before to study the dimerization in the S = 1/2 Heisenberg 
model on a square lattice with an interaction between first and second neighbors |T7 |, and 



also the dimerization transition in a S — 1 chain fl9 |. We believe that this approach has 



some advantages over the mean-field Schwinger-boson theory. For example, it correctly 
reproduces the fact that at the critical point and in the disordered phase, the magnon 
excitation spectrum is three-fold degenerate. 

We now substitute (|) and (|7|) into the Hamiltonian. To leading order in A, the interaction 
between bosons can be neglected, and diagonalizing the quadratic form in bosons, we obtain 
a three-fold degenerate excitation spectrum with the dispersion 



*k = y/Al - Bl (8) 

where A k = J 2 + 2JjV fc , and B k = 2J*u k , and Jj* = J\/X. For sufficiently large J 2 , the 
excitation energy is real (which indicates a stability), and there is a finite gap in the spectrum 



whose minimum is at k = tt. This gap vanishes at J 2 = J 2 = 4 J*. Below this point, the 
excitations near k = tt are purely imaginary which signals an instability and implies a need 
for a change of the ground state. 

To obtain a better estimate for the critical value of J2, we included anharmonic terms 
into consideration, computed the self-energy terms by usual means, and obtained to order 
0(A 2 ): 

J c 2 r = U{ (\ - 0.665A + A 2 logl/A + 0(A 2 )) (9) 

We see that the first-order correction shifts the transition towards smaller J 2- If we had 
restricted the calculation to include only this term, we would obtain J| r in a range between 
1.34Ji and 2.4Ji, depending on whether we leave the correction in the numerator or put 
it into the denominator. The second-order correction is positive and partly compensates 
the downshift renormalization due to the first-order term. Unfortunately, the second-order 



correction is logarithmically divergent at the transition [ 29"| , and we cannot obtain the precise 
value of J 2 to order A 2 . We therefore can only argue that the actual value of J| r is in between 
our zero-order and first-order results. A somewhat better estimate of Jf r can be obtained 
approaching the transition from the ordered phase, and will be discussed in appendix A. 
Notice, however, that the first-order estimate of J% r is already closer to the numerical result 
than ~ 4.3Ji which was obtained in a self-consistent spin-wave approach. 

We also computed the spin- wave velocity at the critical point. To order 0(A), we obtained 

c sw = 2 J* (1 - 0.256 A) (10) 

For the physical case of A = 1, this gives c sw between 1.49Ji and 1.59Ji again depending on 
whether we keep the correction in the numerator or put it into the denominator. The second 
order correction to the spin- wave velocity is again positive and partly compensates the 0(A) 
contribution, but it is again of the form A 2 log 1/ A which prevents us from obtaining the 
precise value of c sw to order 0(A 2 ). Alternatively however, we can reexpress c sw in terms of 
the critical value of J2. Doing this, we find that to order A, c sw = 0.5Jf r (l + 0.409A + ...). 
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For A = 1, this yields c sw = 0.705J| r . Using then the numerical result J| r = 2.55Ji, we 
obtain c sw ~ I.8OJ1 which is consistent with c sw = 1.7 — 1.8 Ji, extracted from the Monte- 
Carlo data. In any case, the velocity we found is smaller than that obtained in the spin-wave 
theory. 



We now consider the case J 2 < J% r when the system possesses a Neel order. We as- 
sume that the sublattice magnetization, N , is directed along the z-axis. In our approach, 
a nonzero iVo = Nq implies that there is a single-particle condensate of the c-quanta 
with momentum n = (7r,7r): < >= a. In a mean-field approximation, we then have 



N = A -1 y/3(l — (3), where (3 = Xa 2 . Introducing the condensate into the Hamiltonian and 
evaluating the ground state energy, E , in the mean-field approximation (i.e., to leading 
order in A, but keeping (3 fixed), we find 



Minimizing the energy, we obtain (3 — (3 — (47^ — J 2 )/8J^. For J 2 = 0, we have (3 = 1/2, 
and hence N = 1/(2A) as it should be. We then performed the standard computations for 
a Bose gas with a condensate and obtained the quasiparticle spectrum. It now contains two 
different branches of quasiparticle excitations. The excitation spectrum for fluctuations in 
the direction perpendicular to the condensate (i.e., for a— and b— type bosons) is doubly 
degenerate. For these excitations, we obtained to leading order in A 



We see that the transverse fluctuations are gapless as they indeed should be. For the spin- 
wave velocity near k = n we have 



III. ORDERED PHASE 




\E = J 2 [3 - 4 J* 13(1 - (3) 



(11) 




(12) 



c. 



2 J* (1 - /3) 



1/2 



(13) 
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Observe that for J2 = we recover the mean-field dispersion for the Heisenberg antiferro- 
magnet: e± = 2J*(1 - z^) 1/2 . 

For the dispersion relation of the fluctuations along the direction of the condensate (i.e., 
for c-type bosons), we found 



We see that the longitudinal fluctuations in the ordered phase have a finite gap at the 



longitudinal mode becomes dispersionless: £\\(k) = 4 J*. However, we do not know whether 
this result survives beyond the leading order in A. The actual dispersion for a c-boson 
may also contain some finite imaginary part (due to higher-order terms in A) which can be 
substantial at small J 2- 

The computations which lead to eq. (|T4|) require some care. The important point is 
that since a ~ A -1 / 2 , there is a cancellation of the overall factor A n in the n-th term in the 
expansion over density in U, and all terms in the series are in fact relevant. In practice, 
this implies that evaluating the contribution to the longitudinal dispersion from L Z L Z , one 
has to examine each term in the series, put all c— bosons except for two into a condensate, 
compute the numerical combinatoric factor, and explicitly sum the resulting series. 

We then used the results for the quasiparticle spectra and computed the sublattice mag- 
netization and the uniform spin susceptibility beyond the mean-field level, to order 0(A). 
The computations and the procedure of extending the first-order results to A = 1 are dis- 
cussed at some length in the appendix. The results are presented in Fig. and Fig. |. For 
comparison, in Fig. |2], we also plotted the self-consistent spin-wave result for the magneti- 
zation. It is essential that at J 2 = 0, both our results are exactly the same as obtained in 
the first-order 1/S expansion. In other words, for a single layer antiferromagnet, there are 
no independent contributions from longitudinal fluctuations. This result provides a qualita- 
tive explanation of why the 1/S" expansion works so well for a single layer antiferromagnet. 
Indeed, in our approach, we treat longitudinal fluctuation as a separate bosonic mode. At 




(14) 



antiferromagnetic momentum, e\\(n) = 8Jj* [(3 (1 — j3)\ 



1//2 . Also observe that at J 2 = 0, the 
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the same time, in the 1/S expansion, the longitudinal mode appears as a pole in the two- 
particle Green function. To obtain this pole, one has to sum an infinite number of the 1/5 
terms. Then, roughly speaking, the contribution from the longitudinal mode represents the 
contributions from high-order terms in the 1/S* expansion. The absence of the longitudi- 
nal correction in our effective "spin- wave theory" therefore implies that the series in 1/S 
converges rapidly, and the dominant contribution comes from the first-order term. 

We emphasize however that the longitudinal fluctuations can be neglected only for 
J2/J1 <S 1- As J2 increases, the deviation of our result for N from the self-consistent 
spin-wave result becomes more and more substantial as seen in Fig. [|. Near the disordering 
transition, longitudinal and transverse fluctuations have nearly equal strength, and the ac- 
tual behavior of sublattice magnetization and uniform susceptibility differs in an essential 
way from the prediction based on the spin-wave theory. 

Notice that in some range of small J2, both the sublattice magnetization and the uniform 
susceptibility are larger than for a single layer, i.e., the system first becomes more "classical", 
and only then, at larger J 2 , do quantum fluctuations push the system towards the disordering 
transition. The region of more "classical" behavior at intermediate J2 has been observed 
in the mean-field Schwinger-boson approach [II]; it is also present in the self-consistent 
spin- wave analysis (see Fig. 

Near the transition point, we obtained 

^0 = ^ v^> XL = A{Z x jU x ) (P/X 2 ) 1/{1+v \ (15) 

where Z N = 1 — 0.163A, Z x = 1 + 0.255A, and 77 ~ 0.03 is the critical exponent for spin 
correlations at criticality. The factor A cannot be obtained within the present approach 
because of the divergence of the Gaussian corrections near the transition point in 2 + 1 
dimensions. Our estimates in the appendix place A to be roughly equal to 2. The ratio 
iVg/[27r(p, ! ) 1 /( 1+,? )] is an overall factor for the dynamical spin susceptibility. Using flI5|) and 
the result for the spin- wave velocity at the transition point, we obtain /[2ir(p s ) 1+v } = 
B/jl +v , where B = (1 — 0.06A)/(27rA 1+r? ). Three different numerical estimates of B all 

11 



yield B = 0.063 [14]]. This is roughly consistent with our estimate B = 0.149/A, though we 
only approximately know that A ~ 2. 

We also computed the quasiparticle dispersion to order 0(A) near the transition, and 
explicitly obtained the Goldstone mode in the transverse channel. These calculations were 
performed only to leading order in /3, when one can neglect cubic terms. For a general 
/3, the Goldstone modes arise as a result of cancellations between the the second-order 
contributions from the cubic terms and the first-order contributions from the quartic terms. 



A similar situation is known to exist in frustrated spin systems ||22j| . We did not perform 
explicit calculations of the spin-wave spectrum at arbitrary (3 and therefore cannot make a 
definite prediction about how longitudinal fluctuations influence the spin-wave velocity at 
small J 2- However, given the good agreement between our result and the spin- wave result 
for the susceptibility in a single layer antiferromagnet, and the consistency between the 
spin-wave result for the spin stiffness, p s = cj. w x±, and the numerical data p0| , we expect 
the corrections due to longitudinal fluctuations to be zero or at least small at vanishing 
J 2 . However, near the transition point, we have already shown that the corrections to the 
spin-wave velocity cannot be reduced to only those due to transverse fluctuations. Thus 
the spin- wave result for c sw , which neglects longitudinal contributions, is most probably 
not quite accurate even though the velocity remains finite at the transition point, and the 
0(1/ S 2 ) correction to c sw is much smaller than the 0(1/ S) correction (see Sec. |I A|) . In other 
words, we argue that near the transition, the series of 1/S terms is not rapidly convergent 
even if the first few terms in the series seem to indicate the contrary. 



IV. UNIFORM SUSCEPTIBILITY AT THE CRITICAL POINT 

In a single- layer Heisenberg antiferromagnet, the linear temperature dependence of the 
uniform susceptibility associated with quantum-critical behavior has been observed in the 
temperature range 0.35Ji < T < 0.6J\. At lower temperatures, there is a crossover to 
another linear behavior associated with the renormalized-classical regime (which, however, 
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has not yet been observed), while at higher temperatures, Xu flattens and has a broad 
maximum at T ~ J\ fl5| , r24}j . How far the linear dependence extends at high T depends 
on the lattice-dependent corrections to scaling. The Monte-Carlo results for a two-layer 
antiferromagnet at the critical J 2 have shown that the linearity extends to sufficiently high 
temperatures, T ~ J 1; i.e., the corrections to scaling at J 2 = J 2 are smaller than those 
of a single-layer antiferromagnet. Below we will compute these corrections perturbatively. 
But first we consider the sigma-model description of a two-layer system, from which one can 
obtain the leading, universal, linear in T term in the uniform susceptibility. 



A. Sigma-model analysis 

A simple way to obtain a sigma-model description of a spin-S" quantum antiferromagnet, 
which we will follow, was suggested by Affleck ||25|| . In application to our system, one has 
to double a unit cell in each of the two layers and introduce n a ^ = (S a> i — S aji+ x)/2S, l a> i = 
(S a ,i + S a) i + i)/2S. At large S, n becomes a classical unit field with commuting components, 
while the commutation relations between n and I are the same as for a vector and a gen- 
erator of rotations. Introducing n a and l a into the Heisenberg Hamiltonian and making a 
transformation from the Hamiltonian to the corresponding action which contains only the 
derivatives of n a , we obtain the action of two interacting 0(3) sigma-models. In terms of ft 
and I, the interaction term has the form S 2 (nin 2 — hh), where 5 2 oc J 2 . The generator of 
rotations itself contains a derivative of n : I ~ n x and the I1I2 term thus only leads to a 
velocity renormalization. Neglecting this term, and also introducing the magnetic field into 
the action for susceptibility calculations, we obtain 

\ 2" 



(Vni) 2 + (Vn 2 ) 2 + S 2 n!n 2 + J- - %H x n\ + J- - iH x n 2 



(16) 



where g is a coupling constant which depends on the ratio J2/J1, and H is measured in units 
of giis/h- Introducing <7 li2 = (n\ ± n 2 )/v / 2, we can rewrite the sigma-model action as 
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1 



(W 1 ) 2 + (W 2 ) 2 + HV 1 2 + 4 



( 



) 



2 




) 



2 



s 



dr 



iH x B\ 



+ 



iH x <j 2 



2^ 



s w 



(17) 



The constraints on the cx-fields are <7i<72 = 0, a\ + cx| = 2. 

The evaluation of the susceptibility at the mean-field (N = oo) level is straightforward. 
Using the results of fU, we obtain Xu = (Xi + X2V2, where x« is a susceptibility per spzn, 
and xi,2 are the mean-field susceptibilities for the two sigma-fields 



where m\ = a/S 2 + m 2 , and rri2 = m, where m is the mass obtained from the second 
constraint equation. At g = g c = 8n(A + \/A 2 + S 2 — |S|) _1 where A ~ J is the upper 
cutoff, we have m = OT + 0(T 2 ), where @ 9 = 21og[( v / 5 + l)/2]. At low T < S, Xi is 
exponentially small in T and can be neglected compared to X2- It is not difficult to show that 
the contributions related to the fluctuations of 0\ are exponentially small and persist even 
beyond the mean-field level. As a result, the universal term in the uniform susceptibility is 
solely due to <t 2 , and Xu is precisely half of that in a single- layer model. 



The sigma-model approach gives us the leading, universal, temperature dependence of 
the uniform susceptibility. Now we compute the leading nonuniversal correction to Xu- We 
will again use a microscopic approach based on a transformation to bosons. However, this 
approach clearly has to be modified compared to what we did before at T = because 
the quasiparticle densities (both normal and anomalous) diverge at finite temperature, and 
the expansion in A is no longer valid. For this reason, we will perform a self-consistent, 
mean-field calculation of the susceptibility: we first assume that anharmonic contributions 
to the quasiparticle spectrum produce a T— dependent gap which eliminates divergencies of 
quasiparticle densities at the transition point, then we evaluate the quasiparticle densities 




(18) 



B. Computation of the subleading term in Xu(T) 
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with the renormalized spectrum, and solve the self-consistent equations for the gap. In 
principle, one can perform these calculations using the same transformation to bosons as 



before. This procedure is then equivalent to self-consistent "1/5"' calculations in 2D |27 
However, we found it more convenient to use a similar but slightly different form of the 
transformation to bosons, introduced by Bhatt and Sachdev ]2?|. In their approach, one 
introduces an extra bosonic field instead of a U— term in ([?]): 

LI = -(4 Si + sfc t ); Lf = V2(ats t + sfh); LJ = V2(bf Si + s+a;). (19) 

The expressions for M are the same as before. The commutation algebra for spins is again 
satisfied, while the constraint on the length of the spin now reduces to afdi + bfbi + cfci + 
sfsi = 1. The advantage of this transformation is that one no longer needs to assume that 
the density of excitations is small. However, we did not use this transformation for our T = 
calculations above because we found it difficult to perform a systematic expansion about the 
mean-field solution. However, the mean-field calculation is straightforward: one has to put 
the s— field into a condensate (< s >= so), neglect fluctuations of s, and reduce the on-site 
constraint to a constraint imposed on average quantities. We first list the T = results which 
are similar (but not identical) to the results we obtained to the zeroth order in A. In the 
disordered phase, we indeed again find the three-fold degenerate quasiparticle spectrum with 



e k = \jA\ — L>1, where A k = J 2 + Uxs\v k , B k = 2J\s\v k , and the self-consistent equation 
for so follows from the constraint on the length of the spin: s\ = 1 — (3/iV) J2k(A k — e k )/2e k . 
At the transition point, we obtained so ~ 0.9. The critical value of J 2 is then J| r = 4JiSq ~ 
3.2J 1; and the T = spin-wave velocity at criticality is c sw = 2J lt s ~ I.8J1. 

We now consider finite temperatures. Assume that the condensate of the s— field has 
a form Sq = (sq)t =0 (1 — m 2 /4), such that at the critical point and near k = e\ = 
c 2 sw {k 2 + m 2 ). Substituting the full expressions for A k and e k into a self-consistency equation 
at finite T, expanding in T and evaluating the lattice sums, we obtain 



e(W| + o(r 2 )) (20) 



T 
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Here O is the same as in the sigma-model calculations, and the second term is a lattice- 
dependent correction which we found to be jj, — —0.061. Furthermore, we have checked 
that the mean-field formula for the uniform susceptibility is given precisely by eq. flI8|) 
with no extra lattice-dependent corrections (we applied a magnetic field, rediagonalized the 
quadratic form in bosons, and computed the magnetization along the field). Substituting 



then the result for the mass m to order T into (Jig), we obtained 

where Q = v^5@/4vr in the mean- field approximation (the 1/N correction extended to a 
physical case of N = 3 reduces this value by about 20% |p3|). We see that the numerical 
factor in the subleading term in the susceptibility is very small, and, e.g., at T = J x , consti- 
tutes only 5% of the mean-field value. Indeed, at T ~ J\, higher-order corrections in T j J\ 
could also be relevant, but the fact that the leading correction to the scaling result is small 
is at least an indication that the universal linear dependence of the uniform susceptibility 
extends to sufficiently high T ~ J\. As we already discussed, this is consistent with the 
Monte-Carlo data O . 



V. CONCLUSIONS 

In this paper, we considered a two-layer Heisenberg antiferromagnet which can either be 
in the Neel-ordered or in the disordered phase at T = depending on the ratio of the intra- 
and interlayer exchange constants. We applied a transformation to bosons which is suitable 
for a singlet configuration of a pair of spins, and considered in a systematic expansion the 
quasiparticle excitations in the disordered phase, and the critical value of the interlayer cou- 
pling. We then extended the approach to the ordered phase by introducing a single-particle 
condensate of one of the Bose fields and computed the mean-field quasiparticle dispersion, 
the sublattice magnetization and the transverse susceptibility at arbitrary J2. We then 
computed one-loop corrections to the sublattice magnetization and the susceptibility, and 
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considered the relative strength of the longitudinal spin fluctuations. We found that the 
contributions of these fluctuations are zero in a single-layer antiferro magnet, but are quite 
substantial near the transition point, where the transverse and the longitudinal fluctuations 
are equally important. The results of our T = calculations are in a reasonable agreement 
with the Monte-Carlo and series expansion data. We also computed the temperature depen- 
dence of the uniform susceptibility at the critical point, and found that the lattice-dependent 
corrections to the universal scaling behavior Xu oc T are small for all T < J\. This is again 
consistent with the Monte-Carlo data which show that the linear behavior of Xu extends to 
sufficiently high temperatures T ~ J\ and flattens only at even higher temperatures. 

It is our pleasure to thank A. Millis, H. Monien, S. Sachdev, A. Sandvik and A. Sokol 
for useful conversations. 



VI. APPENDIX 

In this appendix, we compute the sublattice magnetization, the transverse susceptibility 
and the spin stiffness in the Neel phase to order A. We start with the calculations of the 
magnetization. 



A. Sublattice magnetization 

Our point of departure is the expression for L z , eq. (0), extended to A C 1. In the 
ordered phase, N =< L z > /2 =< c + U > / V\, where the averaging is over the exact 



ground state. The mean-field calculations in the ordered state were presented in Sec. |T|. In 
these calculations, we considered only the condensate piece of the c— field, < c >= a. Here 
we will need both, a and the fluctuating component of c. Substituting Ck = a8k, n + c^, into 
No, expanding in U up to an infinite order, and collecting all terms which contain at most 
one pair product of fluctuating fields, we obtain after some simple combinatorics 

J (3(1 - p) 

N = -^—^ [1 - A (Z t (p) + Z 2 {(3) + Z 3 (J3) + ZM)] (22) 
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where, we recall, (3 = Xa 2 , and 



8(1 -/3) 2 



7(* = .±LlrM 

2{P) 4(1-/3)2 iV^2e||(fe) 



Here 



A\\{k) = J 2 + ((3(4 - 3/3) + u k (l - 2/3) 2 ) , 

B ll(*) = S (/5(2-/3) + ^(l-2/3) 2 ), 

A ± (fe) = J 2 + 2J>* + 4J*/3(1 - u k ), (24) 

•^l = -A/^> an d the dispersions for transverse and longitudinal fluctuations are given by 
(O) and (ITJ]). Observe that near the critical point, iVo oc \ffi. The next step is to express 
(3 in terms of J\ and J 2 . To this end, we compute the ground state energy, E with the 
0(A) corrections which come from noninteracting spin- waves and from the normal ordering 
of c— operators in the expansion of U. Combining the two contributions, we obtain 

XE = J 2 (3 - 4J*/3[(1 -p)-2\(l- (3) Z x {p)\ 

-A YXMk) - e±(k)\ - ± J2[A ]{ (k) - e n (*)] (25) 

k Z k 

Minimization with respect to (3 then yields 

13 = (3 - A (Z 5 (/3 ) + Z 6 (J3q) + Z 7 ((3 ) + Z 8 ((3 )) (26) 

where 



Oj 



1 A 



ZM = £ (4-6/3 + 3/3* + (-3 + 8/3 - 4/3> fc ) (-1 + M| 



ZrW = -^77T^ E (2 - 2/3 + /3 2 + (-3 + 8/3 - 4/?> fc ) ^ ( k ' 



4iV(l-/?) 2 V V - - v - - ' V 2e||(fc) 

Z(B)- P{2 ~ P) (27) 
Zs{f3o) ~ 8(1 -A)) 2 ( } 

Notice that the correction terms Zi, Z 2 , Z 3 , Z 6 , Z 7 and Z 8 are due to fluctuations in the 
direction of the condensate, while the terms Z4 and Z§ come from transverse fluctuations. 



Substituting (f26|) into (P4), we obtain, to order 0(A) 

1 - A Z b /p ^ 1 2 



iV - ^ 



i - A V(i - A>) 



(1 - AZ a ) (28) 



where 



Z a = Zx(A)) + Z 2 (p ) + Z 3 (/3 ) + Z 4 (AO (29) 
Z b = Z 5 (/3 ) + Z 6 (/3 ) + Z 7 (/3 ) + Z 8 (/3 ) 

At J2 = 0, (3q — 1/2, and evaluating the lattice sums, we obtain Nq = (1/2A) — n , where 
n = A^ _1 J]fc((l — V \Y 1 ^ ~~ l)/2 = 0.197 is the density of transverse fluctuations (spin 
waves) p0| . This result is equivalent to the first order spin- wave result, i.e., longitudinal 
fluctuations do not contribute to sublattice magnetization to first order in A. This is a direct 
consequence of the fact that the longitudinal mode is dispersionless at J 2 = 0, and hence 
the c— bosons on adjacent sites do not interact with each other. It is essential, however, that 
the longitudinal fluctuations are small only for J 2 j J\ <C 1. Near the disordering transition, 
longitudinal and transverse fluctuations have nearly equal strength, and the actual behavior 
of magnetization differs in an essential way from the prediction based on the spin-wave 
theory. In this limit, we obtained 

No = -y y/0, (30) 

where Zn = 1 — 0.163A, and the fully renormalized j3 satisfies the equation 

8Ji/3(l - (3/tt) A/y^ + ...) = J 2 cr - J 2 , (31) 

where J 2 = 4J*(1 — 0.665A + . . .) is the same as we obtained approaching the critical point 
from the disordered phase. We have checked that the two analytical expressions for J 2 are 
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indeed also identical. The subleading term in fl3T|) is a Gaussian correction to the sublattice 
magnetization. In the theory of phase transitions, it is usually assumed that the Gaussian 
term is in fact expressed in terms of fully renormalized (3 rather than (3q. The correction 
term then diverges as one approaches J| r as it indeed should in 2+1 dimensions. Due to 
this divergence, the self-consistent approach is valid only at Jf — J 2 > A 2 . In the opposite 
limit J 2 — J 2 <C A 2 , scaling considerations predict that the sublattice magnetization should 
behave as N ~ (Jf r — J2Y , where (3 ~ 0.35. 

The above considerations are also relevant as to how one should extend the perturbative 
result for No at arbitrary /3 to A = 1. We have seen that near J 2 = 0, one should keep 
(3 = (3q in the 0(A) terms. At the same time, it is not difficult to make sure that in order 
to obtain the same J| r on both sides of the transition, one has to perform calculations self- 
consistently, i.e., evaluate the subleading terms in (|26|) with the fully renormalized (3. To 
first order in A, both procedures are indeed equivalent. However, the extension to A = 1 
yields different results in the two cases. The self-consistent solution of fl26|) for A = 1 is 
plotted in Fig. We see that there is a substantial downturn renormalization of (3 in the 
region J| r — J 2 <C A 2 , where the self-consistent solution is in fact invalid. If instead, we 
approximate the critical value of J 2 from the region of intermediate (3 (see Fig. we obtain 
the larger J| r ~ 2.3, which is in better agreement with numerical results. On the other hand, 
the perturbative solution (with (3$ in the subleading terms) gives a correct description of the 
sublattice magnetization at small J 2 , shows no unphysical downturn renormalization near 
the transition, and yields J| r ~ 2.73Ji, which is reasonably close to the numerical result. 
For all these reasons, we plotted the perturbative solution for Nq in Fig. |2|. 

B. Transverse susceptibility 

We will use a direct way to obtain the transverse susceptibility in the ordered phase, 
that is we will apply a homogeneous transverse magnetic field and compute the induced 
magnetization. For definiteness, we will assume that iV is directed along the z— axis, and 
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apply a magnetic field in the x— direction. 

For the calculations of the transverse magnetization, we found it convenient to introduce 
new Bose-operators as linear combinations of the original a— and b— bosons: 

a% + k di- h 

S, = -7T P, = ^T (32) 

In terms of these new operators, the transformation to bosons, extended to large A, is 

Ml = stPi + pfs l , L\ = -\-V 2 (c+Ui + U iCi ), 

Mf = p+ Ci + cfpi, Lf = \- 1/2 (sfU t + Ui8i), 

Mf = -i(stci - c+Si), L\ = -iX-^ipfUi - U iPi ), (33) 

where Ui = (1 - \(sfsi+pfpi + ctd)) 1/2 . 

The advantage of using this new form of the transformation is that a magnetic field 
applied along x, only introduces a condensate of the p— field. As the expectation value of 
M x is obviously site-independent, the c— and p— field condensates should have the same 
momentum, i.e., the condensate of p should also have a momentum tc. 

Let us first discuss the mean-field results. In the mean-field approximation, the transverse 
magnetization per spin is M± = M x /2 =< p > < c >= A -1 (7o/3o) 1//2 , where we have 
introduced 7 = A< p > 2 , and 7 = 70 at the mean-field level. The mean-field ground state 
energy depends on both, (3 Q and 70, and is given by 

XE = J 2 (A) + 7o) - 4Jr/? (l - A)) + 8Jr A)7o - 2H X (70 A)) 1/2 , (34) 

where Eq is the energy per a pair of spins, and, as before, the magnetic field is measured in 
units of gfis/h. Differentiating over 70 and substituting the result into M±, we obtain 

M (7o A)) 1/2 1 (3qH x 

To obtain the susceptibility, we need Mj_ only at vanishing magnetic field. Substituting 
Po = (4 J* - J 2 )/8J* into (H), we find 

1 4J* - J 2 

x± = ui -sir m 
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Observe that at J 2 = 0, we recover the classical spin- wave result x± — 1/8 Ji. For the 



spin-stiffness we obtain using (O) 



Ji 16(Jf) 2 -J, 



Ps 4A 2 16(JJ)2 (37) 
For J2 = 0, we again recover the classical spin- wave result. In the opposite limit, J 2 ~ 4Ji, 
Ps = (Ji/X 2 ) [(AJ*—J 2 )/8Ji\. Finally, for the ratio Nj 2 / p s we have N^/p s = 1 / J x independent 
of J 2 . 

We now obtain the expression for x± to order 0(A). From ([33|) we have Mj_ = X^ 1 y/J3^y + 
AMj_, where AM± = N~ l Y,k^-K < 4 Pfc >• To compute 7 and AM to first order in A, 
we will need the excitation spectra of quasiparticles in the presence of the field. To obtain 
them, we substitute the transformation to bosons into the spin Hamiltonian and restrict 
our calculations to the terms which are quadratic in bosons. The fluctuations of the s— field 
are decoupled from the other two modes, while the fluctuations of the c— and p— fields 
are coupled in the presence of a field. The computation of the quadratic form in the c— 
and p— bosons again requires some care as one needs to carefully examine all terms in the 
expansion of the square root in (p3|), keeping in mind that both 7 and (3 are not small in A. 
Assembling the contributions from all terms in the series, we obtain 

H = E + £ A s {k) s\s k + (4sL fc + s k s. k ) + A p {k)p{p k + ( P yi k + p k p. k ) + 

k ~ ~ 

B (k) 

A c {k) c\c k + (4cL fc + c fc c_ fc ) + C{k)(c\p k +p\c k ) + D(k)(p\cl k +p k c. k ) (38) 

where 

XE = J 2 {(3 + 7) - 4J 1 */?(1 -P)[l- 2XZM] - 2H X ^ + 8 J*/3 7 (l + A/(8(l - /?) 2 ) (39) 
and 

A s (k) = J 2 + 2J{ {2(3 + v k (l -2(3- 2 7 )) , B s {k) = 2J\v k {\ - 2 7 ) 
A v (k) = J 2 + 2J* (2(3 + ^ + ^ fe f 1 - 2(3 - 7 + y^r^ ) , 
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Ac(A;) = j2 + 2 Ji l/ 3 wr + Ff + fe l"T^r - 7 w, 

^ = 2 Jl + Ff + "* - 7 (1-^)2 J J > 

3-2/3-^(2-3/3) 



C{k) = -H x + 2J{ {(3 7 ) 1/2 



/3 



1-^(2-3/3) 



= 2JI 09 7 )Va ^_ ^ (40) 

Diagonalizing then the 2x2 matrix for s— bosons and 4x4 matrix for coupled c— and 
p— bosons by usual means, we obtain three branches of quasiparticle excitations. The dis- 
persion of the s— boson has a gap, e s (7r, tt) = H x . This result is valid for any J2 and is indeed 
the expected result since s-quanta describe the fluctuations of the transverse components of 
M, and the k = tt mode of these fluctuations is a homogeneous precession of the magneti- 
zation around the direction of the field (we recall that both condensates have momentum tt, 
and hence homogeneous (k = 0) modes of composite fields M z and M x correspond to k = tt 
mode of the s— boson). Furthermore, we found after a diagonalization that one of the two 
coupled modes of the c— and p— bosons remain gapless at k = tt, while the other has a gap 
which in the absence of the field is the same as the gap for the c— quanta. The presence of 
a gapless excitation is a direct consequence of the Goldstone theorem. 

For the calculation of the transverse magnetization, we actually need only the ground 
state energy. Collecting the zero-point contributions, which appear after diagonalization, 
we obtain to order H 2 

E tot = Eo- 1 -Y: ElMk) - e^k)} -Y.il ll ^ + ^) + D{l + x p x c )f (41) 
1 i k k e p + 6c 

where E is given by ([39|), i = s,p or c, = (A 2 — B 2 ) 1 ^ 2 , and I 2 = (Ai + ei)/2e$, Xi = 

—Bi/(Ai + £j). Simultaneously, substituting old Bose operators in terms of new ones into 

AMj_, we also obtain 

AMj_ = -Y l l l l ( x p + x c) C,(Xp + Xc) + jP(1 + XpXc) (42) 

k Cp + Cc 



Evaluating 7 oc [3H 2 from dE to tld^ = and using (EH) and (f4"2"D, we obtain the result 



for x± — M±/ H x to order 0(A). The full expression is, however, too cumbersome to be 
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presented here, so we analyze only the limiting cases and plot the result for arbitrary J 2 in 
Fig. (we used the same procedure of extending the result to A = 1 as for the sublattice 
magnetization). Near the critical point, we found 



(7 P) 



1/2 _ pHx 



Jl T - 5\J{Z. 



7 



H„0 /A 2X 1/2 



AM ^=12^fUJ {l+ ° m <43) 
where Z 7 = iV _1 J2k ^fc/v 7 ! + v k = —0.328, J| r is given by ©, and we keep (3 rather than 
(3q in the 0(A) terms. Collecting the two contributions to M±, we obtain 

where Z x = 1 + 0.255A. The subleading term is a Gaussian correction. Its divergence 
again implies that the self-consistent approach only works for (3 > A 2 [29|. In the opposite 



limit, A 2 ^> (3, a self-consistent theory is inapplicable. Scaling considerations 0] predict 
that in this limit x± = A(Z x /4:Jial) (P/X 2 Y^ 1+ri \ where 77 ~ 0.03 is the critical exponent 
for spin correlations at criticality, and A is a constant whose value cannot be obtained in 
the present approach. At A 2 ^> /3, we also have p s = x± c l w = AJiZ p (/5/A 2 ) 1 ^ 14 "^, where 
Z p = l- 0.257A. 



In the opposite limit, J 2 = 0, we find 



XL = ^ («) 

where 

^ = 1-Al E (1 jf |)1/2 =1-0.551A (46) 

is the contribution from s— and p— bosons, which is exactly the same as in the first-order spin- 
wave theory. In other words, longitudinal fluctuations do not contribute to the susceptibility 
of a single-layer antiferromagnet. This is consistent with the fact that the first-order 1/5* 
result for x± agrees well with the numerical data [3~0|. However, the longitudinal fluctuations 



are again small only for J%j J\ <C 1. As J 2 increases, our expression for \u deviates from 
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the spin- wave result, and eventually turns to zero much earlier than in the self-consistent 
spin-wave theory. 
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FIGURES 

FIG. 1. The system under consideration is a two layer antiferromagnet with intralayer exchange 
coupling Ji and interlayer exchange coupling J 2- 

FIG. 2. Sublattice magnetization as a function of J%l J\. Points - the self-consistent spin- wave 
result; solid line - the result of our present calculations which take longitudinal spin fluctuations 
into account. The critical value of interlayer exchange is J| r = 2.73Ji (see appendix A). 

FIG. 3. Transverse susceptibility in the ordered phase as a function of J^j ' J\. The critical value 
of J2 is the same as in Fig.||. 

FIG. 4. The solution of the self-consistent equation for the fully renormalized value of the 
single-particle condensate, j3 =< c n > 2 . Points are the results of the self-consistent calculations 
extended to the physical case of A = 1. The downturn renormalization at small is due to 
divergent Gaussian fluctuations, and is probably unphysical. The solid line is the extrapolation of 
the self-consistent formula at intermediate (3 to (3 = 0. 
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